In 2012, we placed 96 dataloggers recording temperature and humidity across Table Mt. and Silvermine. The loggers were stratified across gradients of elevation, slope, exposure and hillslope position (hilltop, valleybottom), and distance from both False Bay and the Atlantic. In previous analyses, the data has been analyzed to extract diurnal Tmin, Tmax, VPmax (maximum daily vapor pressure), and RHsat.hrs (hours per day with relative humidity above 95%).

All of these analyses and the summary data are available at: https://github.com/dackerly/table-mt-microclimate-2012.git

General summary statistics about the sites and daily weather values can be found in ‘manuscript-summary-stats.r’.

This file explores the results of linear models that analyze spatial variation in daily weather variables. The analyses are run in the accompanying file: TableMt_2012_daily.Rmd, which saves the results to file, to be read in here.

Load files

rm(list=ls())
meta <- read.csv('data/csv_masters/location_meta.csv',as.is=T)
names(meta)
##  [1] "domain"               "siteID"               "logger"              
##  [4] "UTM.east"             "UTM.north"            "LON_DD"              
##  [7] "LAT_DD"               "elevation"            "slope_deg"           
## [10] "aspect_deg"           "magn.true.aspect"     "use4Analyses"        
## [13] "use4ClimateSummaries" "route"                "CollSeq"             
## [16] "series"               "deployed.as"          "deployXssn"          
## [19] "start.date"           "start.hour"           "end.date"            
## [22] "end.hour"             "log.max.d"            "VegSurveyPlot"       
## [25] "VegSurveyDate"        "distOnTransect"       "Orig.siteID"
topo10 <- read.csv('data/csv_masters/topo10.csv',as.is=T)
names(topo10)
##  [1] "siteID"    "domain"    "elevation" "slope"     "aspect"   
##  [6] "d2at"      "d2fb"      "d2cs"      "rad080"    "rad172"   
## [11] "rad355"    "thl0"      "thl315"    "thl337"    "plow050"  
## [16] "plow125"   "plow250"   "plow500"   "tpi050"    "tpi125"   
## [21] "tpi250"    "tpi500"
dlySummary <- read.csv('data/csv_outfiles/dlySummary.csv',as.is=T)
head(dlySummary)
##   doy year month day       date     Tmin     Tmax    Tmean   dtRange
## 1   1 2012    NA  NA 2012-01-01 11.39483 21.62492 16.50988 10.230091
## 2   2 2012    NA  NA 2012-01-02 14.97574 19.01598 16.99586  4.040239
## 3   3 2012    NA  NA 2012-01-03 14.75741 23.10118 18.92930  8.343773
## 4   4 2012    NA  NA 2012-01-04 10.78530 23.34831 17.06680 12.563011
## 5   5 2012    NA  NA 2012-01-05 13.78350 25.12926 19.45638 11.345761
## 6   6 2012    NA  NA 2012-01-06 14.36026 27.59257 20.97641 13.232307
##        T06      T14   Tm1406   dt1406    RHmin    RHmax    VPmax RHsat.hrs
## 1 12.01950 20.73315 16.37632 8.713648 83.64332 92.68127 1.920619  5.545455
## 2 15.26168 18.31991 16.79080 3.058227 52.31199 99.71526 1.937674 19.000000
## 3 15.04712 22.19107 18.61910 7.143943 52.63628 98.78947 1.702166 15.679924
## 4 12.20999 21.74176 16.97587 9.531773 39.66210 94.06428 1.792897  4.657197
## 5 14.83984 24.38744 19.61364 9.547602 28.20923 95.26178 1.872311  7.526515
## 6 16.17932 24.65658 20.41795 8.477261 51.95174 85.88167 2.043691  1.242424
##   TminSr TmaxSr   RHsatSr   VPmaxSr sommax sommin Kps.mean_hPa Kppt_mm
## 1 11.897 15.774 17.666667 0.8270558     16     15        997.0     0.6
## 2  5.775 13.309 14.666667 0.6197326     15     14        995.5    11.6
## 3  7.611 14.073 24.000000 0.6384510      7      7        997.2     0.2
## 4 13.317 12.957 13.000000 0.6167037     11     11        998.2     0.0
## 5  7.701  9.869 14.333333 0.5638149     11      7        999.6     0.0
## 6  9.235  9.813  7.833333 0.5298096     11      8        999.4     0.0
##   Kwspd.mean_m.s Kwspd.max_m.s
## 1       2.991667           5.6
## 2       2.504167           5.4
## 3       2.279167           4.7
## 4       2.304167           3.1
## 5       1.975000           2.7
## 6       2.170833           4.8
vars <- c('Tmin','Tmax','T06','T14','RHsat.hrs','VPmax')

Load results

In the previous script, I ran linear models for each day of the year to explain spatial variation in daily weather data, for six variables: Tmin, Tmax, T06, T14, RHsat.hrs, and VPmax.

Models were run using elevation only (just the lapse rate effect). Then a set of 5-term models were run that swapped alternative variables for three different parameters, using all possible combinations, for a total of 168 models. The five parameters were:

1: elevation 2: solar radiation, with 7 alternative terms: (‘dsol’,‘rad080’,‘rad172’,‘rad355’,‘thl315’,‘thl337’,‘thl0’) 3: hillslope position, with 8 alternative terms: (‘tpi050’,‘tpi125’,‘tpi250’,‘tpi500’,‘plow050’,‘plow125’,‘plow250’,‘plow500’) 4: regional location in terms of distance to False Bay or the Atlantic Ocean: (‘d2fb’,‘d2at’,‘d2cs’) 5: ‘slope’

The results are stored in 9 lists, each one containing four items with a matrix for each of the four weather variables RWQelev: r-squared for elevation only model (vector, 366) MSEelev: mean-square error for elevation only model (vector, 366) RSQm5: r-squared for all 5-term models (matrix, 366 x 168) MSEm5: mean-square error for all 5-term models (matrix, 366 x 168) RSQm5x: highest r-squared for 5-term models (vector, 366) MSEm5x: MSE for model with highest R-square, of 5-term models (vector, 366) BGm5x: model number for best fit model (vector, 366) modterms: matrix of terms for each model (matrix, 168 x 5) modnums: matrix of term ids for each model (matrix, 168 x 5)

allRes <- readRDS('data/Rdata/M5out.Rdata')
RSQelev <- allRes[[1]]
MSEelev <- allRes[[2]]
RSQm5 <- allRes[[3]]
MSEm5 <- allRes[[4]]
RSQm5x <- allRes[[5]]
MSEm5x <- allRes[[6]]
BFm5x <- allRes[[7]]
dAICx <- allRes[[8]]
rSlpx <- allRes[[9]]
modterms <- allRes[[10]]
modnums <- allRes[[11]]

vars <- c('Tmin','Tmax','T06','T14','RHsat.hrs','VPmax')

Here are the first 6 models, to show what’s in the modterms matrix. The 168 rows loop through the 3 variable terms (columns 2-4) in the order shown above.

dim(modterms)
## [1] 168   5
head(modterms)
##      [,1]        [,2]   [,3]     [,4]   [,5]   
## [1,] "elevation" "dsol" "tpi050" "d2fb" "slope"
## [2,] "elevation" "dsol" "tpi050" "d2at" "slope"
## [3,] "elevation" "dsol" "tpi050" "d2cs" "slope"
## [4,] "elevation" "dsol" "tpi125" "d2fb" "slope"
## [5,] "elevation" "dsol" "tpi125" "d2at" "slope"
## [6,] "elevation" "dsol" "tpi125" "d2cs" "slope"

Results

Now, let’s look at the results! First, how do Rsq for the elevation and the best of the 5-term models vary through the year, for each variable? Interestingly, for all four variables results fluctuate a lot day to day, especially based on elevation alone. Using 5 variables, r-squareds are mostly between 40 and 90%, with lower values for: Tmin in spring and fall, Tmax in summer, RHsat.hrs in winter, and VPmax in summer.

op=par(mfrow=c(3,2),mar=c(5,5,2,1))
i=1
for (i in 1:length(vars)) {
    plot(1:366,RSQelev[[i]],type='l',col='red',ylim=c(0,1),xlab='DOY',ylab='r-squared',main=paste('R2 of 5-term vs. elevation only:',vars[i]))
    points(1:366,RSQm5x[[i]],type='l',lwd=2)
}

par(op)

Plots of Rsquared for elevation only vs. the five term models are more useful to see how much the other terms can add, especially when elevation alone is a weak predictor. But you can also see days on which the r2 for both models is low - none of these predictors are sufficient to capture the observed variation. We return to this below.

op=par(mfrow=c(3,2),mar=c(5,5,2,1))
i=1
for (i in 1:length(vars)) {
    plot(RSQm5x[[i]]~RSQelev[[i]],ylim=c(0,1),xlim=c(0,1),main=vars[i],xlab='Elevation model r2',ylab='Full model r2')
    abline(0,1)
}

par(op)

Interestingly, there are only weak correlations between model fit across variables on a given day. The strongest correlation is for the 5-term models for Tmax and Tm14 (var 2 vs var 4 in second plot below):

pairs(cbind(RSQelev[[1]],RSQelev[[2]],RSQelev[[3]],RSQelev[[4]],RSQelev[[5]],RSQelev[[6]]))

pairs(cbind(RSQm5x[[1]],RSQm5x[[2]],RSQm5x[[3]],RSQm5x[[4]],RSQm5x[[5]],RSQm5x[[6]]))

My main goal has been trying to make sense of variation in the daily r-squared variables, and next step would be to sort out the meaning of which terms contribute to the best model. Presumably many of the r-squared values are effectively indistinguishable, but that would require quite a lot of examination.

Here are plots of the max r-squared for each variable versus the daily mean for that variable.

Tmin - nothing doing!

plot(dlySummary$Tmin,RSQm5x[[1]],ylab='Max r-squared, Tmin')

Tmax - r-squared are higher on cooler days, i.e. topography is more predictive of spatial variation in Tmax

plot(dlySummary$Tmax,RSQm5x[[2]],ylab='Max r-squared, Tmax')

Tmax variation is also much better explained on days with a low Tmin-Tmax difference, which would presumably be humid and/or cloudy:

plot(dlySummary$Tmax-dlySummary$Tmin,RSQm5x[[2]],ylab='Max r-squared, Tmax',xlab='Tmax-Tmin')

That’s interesting, as I think of the coolest days as often being clear with a potentially high Tmin-Tmax difference. Just a quick check on how those two are related:

plot(dlySummary$Tmax,dlySummary$Tmax-dlySummary$Tmin,xlab='Tmax',ylab='Tmax-Tmin')

plot(dlySummary$Tmin,dlySummary$Tmax,xlab='Tmin',ylab='Tmax')
abline(0,1)

T06 - a little more than Tmin - strongest r-squared at colder temperatures, but not all cooler days

plot(dlySummary$T06,RSQm5x[[3]],ylab='Max r-squared, Tmin')

T14 - even stronger than Tmax - higher r-squared on days with cooler afternoon temperatures

plot(dlySummary$T14,RSQm5x[[4]],ylab='Max r-squared, Tmin')

RHsat.hrs - r-squared are higher on days with intermediate number of hours of fog - perhaps these days have the most variation among stations to explain, since clear (=0) and all foggy (=24), there’s not much variation to work with

plot(dlySummary$RHsat.hrs,RSQm5x[[5]],ylab='Max r-squared, RHsat.hrs')

VPmax - not much going on

plot(dlySummary$VPmax,RSQm5x[[6]],ylab='Max r-squared, VPmax')

Another possibility is that the r-squareds are higher on days with a greater range of values, just because there’s more variation to explain. These plots show r-squared for the elevation model and then the 5-term model, vs. the range from min to max across all 90 stations for a given day. I haven’t calculated the station range for T06 and T14 so those are not shown here.

op=par(mfrow=c(2,2),mar=c(5,5,1,1))
plot(dlySummary$TminSr,RSQelev[[1]],ylab='Elev r-squared, Tmin',xlab='Range of Tmin across sites')
plot(dlySummary$TmaxSr,RSQelev[[2]],ylab='Elev r-squared, Tmax',xlab='Range of Tmax across sites')
plot(dlySummary$RHsatSr,RSQelev[[5]],ylab='Elev r-squared, RHsat.hrs',xlab='Range of RHsat across sites')
plot(dlySummary$VPmaxSr,RSQelev[[6]],ylab='Elev r-squared, VPmax',xlab='Range of VPmax across sites')

par(op)

op=par(mfrow=c(2,2),mar=c(5,5,1,1))
plot(dlySummary$TminSr,RSQm5x[[1]],ylab='M5 r-squared, Tmin',xlab='Range of Tmin across sites')
plot(dlySummary$TmaxSr,RSQm5x[[2]],ylab='M5 r-squared, Tmax',xlab='Range of Tmax across sites')
plot(dlySummary$RHsatSr,RSQm5x[[5]],ylab='M5 r-squared, RHsat.hrs',xlab='Range of RHsat across sites')
plot(dlySummary$VPmaxSr,RSQm5x[[6]],ylab='M5 r-squared, VPmax',xlab='Range of VPmax across sites')

par(op)

Well, that’s interesting. It’s the opposite patterns of what was expected, at least for Tmin. The spatial models are best at predicting variation when the range is low. Note that the monitoring sites span about 800 m elevational range, so a perfect environmental lapse rate would predict about 5°C variation. It’s at that low range around 5° that the models have the highest R-squared, suggesting that higher ranges occur when other factors overwhelm the lapse rate effect.

There’s also not much obvious variation among SOM classes from Jasper’s regional weather pattern analysis:

plot(dlySummary$sommin,RSQm5x[[1]],ylab='Max r-squared, Tmin',xlab='Tmin regional SOM')

plot(dlySummary$sommax,RSQm5x[[2]],ylab='Max r-squared, Tmax',xlab='Tmax regional SOM')

plot(dlySummary$sommin,RSQm5x[[3]],ylab='Max r-squared, T06',xlab='Tmin regional SOM')

plot(dlySummary$sommax,RSQm5x[[4]],ylab='Max r-squared, T14',xlab='Tmax regional SOM')

The model fitting process chose the best term from each set, but perhaps more important is the contribution of each term to the 5 term model, through the year. This was captured with deltaAIC values from a drop1 function on the best-fit model.

Here are the seasonal trends in importance for the five terms, for Tmin, followed by the trends in the actual slopes of the regression terms for each parameter

op=par(mfrow=c(3,2),mar=c(5,5,1,1))
plot(1:366,dAICx[[1]][,1],type='l',ylab='Tmin-elev dAIC',xlab='day of year')
plot(1:366,dAICx[[1]][,2],type='l',ylab='Tmin-solrad dAIC',xlab='day of year')
plot(1:366,dAICx[[1]][,3],type='l',ylab='Tmin-hill dAIC',xlab='day of year')
plot(1:366,dAICx[[1]][,4],type='l',ylab='Tmin-oc.dist dAIC',xlab='day of year')
plot(1:366,dAICx[[1]][,5],type='l',ylab='Tmin-slope dAIC',xlab='day of year')
par(op)

#pairs(dAICx[[1]])

op=par(mfrow=c(3,2),mar=c(5,5,1,1))
plot(1:366,rSlpx[[1]][,1],type='l',ylab='Tmin-elev slope',xlab='day of year')
abline(h=c(0,-0.006),lty=c(2,3))
plot(1:366,rSlpx[[1]][,2],type='l',ylab='Tmin-solrad slope',xlab='day of year')
abline(h=0,lty=2)
plot(1:366,rSlpx[[1]][,3],type='l',ylab='Tmin-hill slope',xlab='day of year')
abline(h=0,lty=2)
plot(1:366,rSlpx[[1]][,4],type='l',ylab='Tmin-oc.dist slope',xlab='day of year')
abline(h=0,lty=2)
plot(1:366,rSlpx[[1]][,5],type='l',ylab='Tmin-slope slope',xlab='day of year')
abline(h=0,lty=2)
par(op)

#pairs(dAICx[[1]])

# 
# op=par(mfrow=c(3,2),mar=c(5,5,1,1))
# plot(rSlpx[[1]][,1],dAICx[[1]][,1])
# plot(rSlpx[[1]][,2],dAICx[[1]][,2])
# plot(rSlpx[[1]][,3],dAICx[[1]][,3])
# plot(rSlpx[[1]][,4],dAICx[[1]][,4])
# plot(rSlpx[[1]][,5],dAICx[[1]][,5])
# par(op)
# pairs(rSlpx[[1]])

For Tmin, it appears that environmental lapse rates fluctuate around expected values of about -6°/1000 m (-0.006 here), with slightly steeper rates in spring. Other terms have expected effects: positive effect of solar radiation (units depends on exactly which term was selected in the model);

Here are the seasonal trends in importance and the slopes for the five terms, for Tmax

op=par(mfrow=c(3,2),mar=c(5,5,1,1))
plot(1:366,dAICx[[2]][,1],type='l',ylab='Tmax-elev dAIC',xlab='day of year')

plot(1:366,dAICx[[2]][,2],type='l',ylab='Tmax-solrad dAIC',xlab='day of year')
plot(1:366,dAICx[[2]][,3],type='l',ylab='Tmax-hill dAIC',xlab='day of year')
plot(1:366,dAICx[[2]][,4],type='l',ylab='Tmax-oc.dist dAIC',xlab='day of year')
plot(1:366,dAICx[[2]][,5],type='l',ylab='Tmax-slope dAIC',xlab='day of year')
par(op)

#pairs(dAICx[[1]])

op=par(mfrow=c(3,2),mar=c(5,5,1,1))
plot(1:366,rSlpx[[2]][,1],type='l',ylab='Tmax-elev slope',xlab='day of year')
abline(h=c(0,-0.006),lty=c(2,3))
plot(1:366,rSlpx[[2]][,2],type='l',ylab='Tmax-solrad slope',xlab='day of year')
abline(h=0,lty=2)
plot(1:366,rSlpx[[2]][,3],type='l',ylab='Tmax-hill slope',xlab='day of year')
abline(h=0,lty=2)
plot(1:366,rSlpx[[2]][,4],type='l',ylab='Tmax-oc.dist slope',xlab='day of year')
abline(h=0,lty=2)
plot(1:366,rSlpx[[2]][,5],type='l',ylab='Tmax-slope slope',xlab='day of year')
abline(h=0,lty=2)
par(op)

#pairs(dAICx[[1]])

Here are the seasonal trends in importance and the slopes for the five terms, for T06

op=par(mfrow=c(3,2),mar=c(5,5,1,1))
plot(1:366,dAICx[[3]][,1],type='l',ylab='T06-elev dAIC',xlab='day of year')

plot(1:366,dAICx[[3]][,2],type='l',ylab='T06-solrad dAIC',xlab='day of year')
plot(1:366,dAICx[[3]][,3],type='l',ylab='T06-hill dAIC',xlab='day of year')
plot(1:366,dAICx[[3]][,4],type='l',ylab='T06-oc.dist dAIC',xlab='day of year')
plot(1:366,dAICx[[3]][,5],type='l',ylab='T06-slope dAIC',xlab='day of year')
par(op)

#pairs(dAICx[[1]])

op=par(mfrow=c(3,2),mar=c(5,5,1,1))
plot(1:366,rSlpx[[3]][,1],type='l',ylab='T06-elev slope',xlab='day of year')
abline(h=c(0,-0.006),lty=c(2,3))
plot(1:366,rSlpx[[3]][,2],type='l',ylab='T06-solrad slope',xlab='day of year')
abline(h=0,lty=2)
plot(1:366,rSlpx[[3]][,3],type='l',ylab='T06-hill slope',xlab='day of year')
abline(h=0,lty=2)
plot(1:366,rSlpx[[3]][,4],type='l',ylab='T06-oc.dist slope',xlab='day of year')
abline(h=0,lty=2)
plot(1:366,rSlpx[[3]][,5],type='l',ylab='T06-slope slope',xlab='day of year')
abline(h=0,lty=2)
par(op)

#pairs(dAICx[[1]])

Here are the seasonal trends in importance and the slopes for the five terms, for T14. Elevation, solar radiation and slope are much stronger in winter. Slopes for ocean distance drop during winter.

op=par(mfrow=c(3,2),mar=c(5,5,1,1))
plot(1:366,dAICx[[4]][,1],type='l',ylab='T14-elev dAIC',xlab='day of year')

plot(1:366,dAICx[[4]][,2],type='l',ylab='T14-solrad dAIC',xlab='day of year')
plot(1:366,dAICx[[4]][,3],type='l',ylab='T14-hill dAIC',xlab='day of year')
plot(1:366,dAICx[[4]][,4],type='l',ylab='T14-oc.dist dAIC',xlab='day of year')
plot(1:366,dAICx[[4]][,5],type='l',ylab='T14-slope dAIC',xlab='day of year')
par(op)

#pairs(dAICx[[1]])

op=par(mfrow=c(3,2),mar=c(5,5,1,1))
plot(1:366,rSlpx[[4]][,1],type='l',ylab='T14-elev slope',xlab='day of year')
abline(h=c(0,-0.006),lty=c(2,3))
plot(1:366,rSlpx[[4]][,2],type='l',ylab='T14-solrad slope',xlab='day of year')
abline(h=0,lty=2)
plot(1:366,rSlpx[[4]][,3],type='l',ylab='T14-hill slope',xlab='day of year')
abline(h=0,lty=2)
plot(1:366,rSlpx[[4]][,4],type='l',ylab='T14-oc.dist slope',xlab='day of year')
abline(h=0,lty=2)
plot(1:366,rSlpx[[4]][,5],type='l',ylab='T14-slope slope',xlab='day of year')
abline(h=0,lty=2)
par(op)

#pairs(dAICx[[1]])

Here are the seasonal trends in importance for the five terms, for RHsat.hrs

op=par(mfrow=c(3,2),mar=c(5,5,1,1))
plot(1:366,dAICx[[5]][,1],type='l',ylab='RHsat-elev dAIC',xlab='day of year')

plot(1:366,dAICx[[5]][,2],type='l',ylab='RHsat-solrad dAIC',xlab='day of year')
plot(1:366,dAICx[[5]][,3],type='l',ylab='RHsat-hill dAIC',xlab='day of year')
plot(1:366,dAICx[[5]][,4],type='l',ylab='RHsat-oc.dist dAIC',xlab='day of year')
plot(1:366,dAICx[[5]][,5],type='l',ylab='RHsat-slope dAIC',xlab='day of year')
par(op)

#pairs(dAICx[[1]])

op=par(mfrow=c(3,2),mar=c(5,5,1,1))
plot(1:366,rSlpx[[5]][,1],type='l',ylab='RHsat-elev slope',xlab='day of year')
abline(h=c(0,-0.006),lty=c(2,3))
plot(1:366,rSlpx[[5]][,2],type='l',ylab='RHsat-solrad slope',xlab='day of year')
abline(h=0,lty=2)
plot(1:366,rSlpx[[5]][,3],type='l',ylab='RHsat-hill slope',xlab='day of year')
abline(h=0,lty=2)
plot(1:366,rSlpx[[5]][,4],type='l',ylab='RHsat-oc.dist slope',xlab='day of year')
abline(h=0,lty=2)
plot(1:366,rSlpx[[5]][,5],type='l',ylab='RHsat-slope slope',xlab='day of year')
abline(h=0,lty=2)
par(op)

#pairs(dAICx[[1]])

Here are the seasonal trends in importance for the five terms, for VPmax

op=par(mfrow=c(3,2),mar=c(5,5,1,1))
plot(1:366,dAICx[[6]][,1],type='l',ylab='VPmax-elev dAIC',xlab='day of year')

plot(1:366,dAICx[[6]][,2],type='l',ylab='VPmax-solrad dAIC',xlab='day of year')
plot(1:366,dAICx[[6]][,3],type='l',ylab='VPmax-hill dAIC',xlab='day of year')
plot(1:366,dAICx[[6]][,4],type='l',ylab='VPmax-oc.dist dAIC',xlab='day of year')
plot(1:366,dAICx[[6]][,5],type='l',ylab='VPmax-slope dAIC',xlab='day of year')
par(op)

#pairs(dAICx[[1]])

op=par(mfrow=c(3,2),mar=c(5,5,1,1))
plot(1:366,rSlpx[[6]][,1],type='l',ylab='VPmax-elev slope',xlab='day of year')
abline(h=c(0,-0.006),lty=c(2,3))
plot(1:366,rSlpx[[6]][,2],type='l',ylab='VPmax-solrad slope',xlab='day of year')
abline(h=0,lty=2)
plot(1:366,rSlpx[[6]][,3],type='l',ylab='VPmax-hill slope',xlab='day of year')
abline(h=0,lty=2)
plot(1:366,rSlpx[[6]][,4],type='l',ylab='VPmax-oc.dist slope',xlab='day of year')
abline(h=0,lty=2)
plot(1:366,rSlpx[[6]][,5],type='l',ylab='VPmax-slope slope',xlab='day of year')
abline(h=0,lty=2)
par(op)

# t1 <- table(BFm5x[[1]])
# modterms[as.numeric(names(t1)[which.max(t1)]),]
# 
# t1 <- table(BFm5x[[2]])
# modterms[as.numeric(names(t1)[which.max(t1)]),]
# 
# t1 <- table(BFm5x[[3]])
# modterms[as.numeric(names(t1)[which.max(t1)]),]
# 
# t1 <- table(BFm5x[[4]])
# modterms[as.numeric(names(t1)[which.max(t1)]),]

# modterms[46,]
# modnums[46,]
# plot(modnums[BFm5x[[2]],3])
# 
# plot(RSQelev[[1]],RSQm5x[[1]])
# plot(RSQelev[[2]],RSQm5x[[2]])
# plot(RSQelev[[3]],RSQm5x[[3]])
# plot(RSQelev[[4]],RSQm5x[[4]])
# 
# plot(MSEelev[[1]],MSEm5x[[1]]);abline(0,1)
# plot(MSEelev[[2]],MSEm5x[[2]]);abline(0,1)
# plot(MSEelev[[3]],MSEm5x[[3]]);abline(0,1)
# plot(MSEelev[[4]],MSEm5x[[4]]);abline(0,1)
# 
# head(dlySummary)
# boxplot(MSEm5x[[2]]~dlySummary$sommax)
# boxplot(MSEm5x[[1]]~dlySummary$sommin)
# 
# plot(RSQm5x[[1]]~dlySummary$Kps.mean_hPa)
# plot(RSQm5x[[2]]~dlySummary$Kps.mean_hPa)
# plot(RSQm5x[[3]]~dlySummary$Kps.mean_hPa)
# plot(RSQm5x[[4]]~dlySummary$Kps.mean_hPa)
# 
# plot(RSQm5x[[1]]~dlySummary$Kwspd.mean_m.s)
# plot(RSQm5x[[2]]~dlySummary$Kwspd.mean_m.s)
# plot(RSQm5x[[3]]~dlySummary$Kwspd.mean_m.s)
# plot(RSQm5x[[4]]~dlySummary$Kwspd.mean_m.s)
# 
# plot(RSQm5x[[1]]~dlySummary$Tmin)
# plot(RSQm5x[[2]]~dlySummary$Tmax)
# plot(RSQm5x[[3]]~dlySummary$RHsat.hrs)
# plot(RSQm5x[[4]]~dlySummary$VPmax)
# 
# plot(RSQm5x[[1]]~I(dlySummary$Tmax - dlySummary$Tmin))
# plot(RSQm5x[[2]]~I(dlySummary$Tmax - dlySummary$Tmin))
# plot(RSQm5x[[3]]~I(dlySummary$Tmax - dlySummary$Tmin))
# plot(RSQm5x[[4]]~I(dlySummary$Tmax - dlySummary$Tmin))
# 
# plot(I(dlySummary$Tmax - dlySummary$Tmin)~dlySummary$Tmin)
# plot(I(dlySummary$Tmax - dlySummary$Tmin)~dlySummary$Tmax)
# plot(I(dlySummary$Tmax - dlySummary$Tmin)~dlySummary$RHsat.hrs)
# plot(I(dlySummary$Tmax - dlySummary$Tmin)~dlySummary$VPmax)